
# 加载所需库------
if (!require("pacman")) install.packages("pacman")
pkgs <- c("dplyr", "ggplot2", "readxl", "tidyr", "lubridate", "janitor", "stringr", "writexl",
          "psych", "magrittr", "sf", "spdep", "tmap", "ggspatial","RColorBrewer")
pacman::p_load(pkgs, character.only = TRUE)

# 读取数据并初步处理
data <- read_excel("./data/processed/cleaned_data9.xlsx") %>%
  filter(year >= 2005) %>% 
  mutate(birth_date = ymd(birth_date), 
         onset_date = ymd(onset_date),
         diag_time = ymd(diag_time))

# 加载人口数据
source("./code/2.3 人口数估发病率.R")

head(data)
is.na(data) %>% sum()
summary(data)
data$city_code %>% unique()
names(data)

# 定义年龄组的水平和顺序
age_levels <- c("0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "≥70")

# 将 age_group 转换为有序因子
data <- data %>%
  mutate(age_group = factor(age_group, levels = age_levels, ordered = TRUE))
names(data)

# 定义一个通用函数来计算分布及占比
calculate_distribution <- function(df, group_cols, count_col_name = "count") {
  # 检查提供的列名是否存在于数据框中
  missing_cols <- setdiff(group_cols, colnames(df))
  if (length(missing_cols) > 0) {
    stop(paste("The following columns do not exist in the dataframe:", paste(missing_cols, collapse = ", ")))
  }
  
  df %>%
    count(across(all_of(group_cols)), name = count_col_name) %>%
    group_by(across(all_of(group_cols[-length(group_cols)]))) %>%
    mutate(total = sum(!!sym(count_col_name), na.rm = TRUE),
           percentage = (!!sym(count_col_name)) / total * 100,
           .groups = 'drop') %>%
    ungroup() %>% 
    select(-.groups)  # 确保移除 .groups 列
}

# 使用字符串直接指定列名
gender_distribution <- calculate_distribution(data, c("sex"))
gender_distribution

# 1. 人群分布特征分析------

## 1.1 性别特征------

### 1.1.1 总体发病情况和分性别发病数及占比------
gender_distribution <- calculate_distribution(data, c("sex"))

print("总体发病情况和分性别发病数及占比:")
print(gender_distribution)
# writexl::write_xlsx(gender_distribution, "./output/表格1.1.1_人群分布_性别.xlsx")

wb <- createWorkbook()  # 新建空工作簿（无任何 Sheet）
if (!"性别" %in% names(wb)) {addWorksheet(wb, "性别")}
writeData(wb, sheet = "性别", x = gender_distribution)

### 1.1.2 不同年龄阶段的性别比------
age_gender_ratio <- data %>%
  group_by(age_group, sex) %>%
  summarise(count = n(), .groups = 'drop') %>%
  group_by(age_group) %>%
  mutate(total = sum(count),
         ratio = count / total)

print("不同年龄阶段的性别比:")
print(age_gender_ratio)
writexl::write_xlsx(age_gender_ratio, "./output/表格1.1.2_人群分布_年龄.xlsx")
if (!"年龄" %in% names(wb)) {addWorksheet(wb, "年龄")}
writeData(wb, sheet = "年龄", x = age_gender_ratio)

### 1.1.3 分年份男女占比情况------
yearly_gender_dist <- calculate_distribution(data, c("year", "sex")) %>%
  mutate(sex = case_when(
    sex == "female" ~ "Female",
    sex == "male" ~ "Male",
    TRUE ~ sex
  ))

print("分年份男女占比情况:")
print(yearly_gender_dist)
writexl::write_xlsx(yearly_gender_dist, "./output/表格1.1.3_人群分布_分年份男女.xlsx")

if (!"逐年男女" %in% names(wb)) {addWorksheet(wb, "逐年男女")}
writeData(wb, sheet = "逐年男女", x = yearly_gender_dist)

# 绘制堆叠柱状图
ggplot(yearly_gender_dist, aes(x = factor(year), y = percentage, fill = sex)) +
  geom_col(width = 0.7) +  # 控制柱子宽度
  labs(x = "Year", 
       y = "Percentage (%)", 
       fill = "Gender",
       title = "Gender Distribution by Year") +
  scale_fill_manual(values = c("Female" = "#FF6B6B", "Male" = "#4ECDC4")) +  # 自定义颜色
  scale_y_continuous(expand = c(0, 0)) +  # 消除y轴空白
  theme_minimal(base_size = 12) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),  # 倾斜x轴标签
    panel.grid.major.x = element_blank(),  # 去除垂直网格线
    legend.position = "top"
  )
ggsave("./figs/图1 分性别发病构成比.png")

### 1.1.4 分月份男女占比情况------
monthly_gender_dist <- calculate_distribution(data, c("month", "sex"))

print("分月份男女占比情况:")
print(monthly_gender_dist)
writexl::write_xlsx(monthly_gender_dist, "./output/表格1.1.4_人群分布_分月份男女.xlsx")

if (!"逐月男女" %in% names(wb)) {addWorksheet(wb, "逐月男女")}
writeData(wb, sheet = "逐月男女", x = monthly_gender_dist)


## 1.2 年龄特征------

### 1.2.1 各年龄组发病情况（发病数）及占比------
age_group_dist <- calculate_distribution(data, c("age_group"))

print("各年龄组发病情况（发病数）及占比:")
print(age_group_dist)
writexl::write_xlsx(age_group_dist, "./output/表格1.2.1_人群分布_各年龄组.xlsx")


if (!"年龄分布" %in% names(wb)) {addWorksheet(wb, "年龄分布")}
writeData(wb, sheet = "年龄分布", x = age_group_dist)


### 1.2.2 分年份各年龄组发病情况（发病数）及占比------
yearly_age_group_dist <- calculate_distribution(data, c("year", "age_group")) %>%
  arrange(year, age_group) %>% filter(!is.na(age_group))
writexl::write_xlsx(yearly_age_group_dist, "./output/表格1.2.2_人群分布_分年份各年龄组.xlsx")

if (!"年龄分布_逐年" %in% names(wb)) {addWorksheet(wb, "年龄分布_逐年")}
writeData(wb, sheet = "年龄分布_逐年", x = yearly_age_group_dist)


# 绘制图1: 分年份各年龄组发病情况及占比
# 查看调色板颜色

RColorBrewer::display.brewer.pal(n = 8, name = "Set2")  # 显示8色Set2方案
p_各年龄组发病情况及占比 <- ggplot(yearly_age_group_dist, aes(x = factor(year), y = percentage, fill = age_group)) +
  geom_bar(stat = "identity") +
  labs(
    # title = "各年份各年龄组发病构成比",
       x = "Year",
       y = "Percentage (%)",
       fill = "Age Group") +
  theme_minimal() +
  scale_fill_brewer(palette = "Set3") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom")#+
  # guides(fill = guide_legend(reverse = TRUE))  # 反转图例顺序，但不影响堆叠顺序
# ggplot(yearly_age_group_dist, aes(x = factor(year), y = percentage, fill = age_group)) +
#   geom_bar(stat = "identity") +
#   labs(
#     title = "Age-specific Incidence Proportions Over Time",
#     x = "Year",
#     y = "Proportion (%)",
#     fill = "Age Group"
#   ) +
#   theme_minimal(base_size = 12) +
#   theme(
#     axis.text.x = element_text(angle = 45, hjust = 1, color = "black"),
#     plot.title = element_text(face = "bold", hjust = 0.5),
#     legend.position = "bottom"
#   ) +
#   scale_fill_brewer(
#     palette = "Set2",  # 使用RColorBrewer调色板
#     # breaks = rev(levels(yearly_age_group_dist$age_group)),  # 反转图例顺序
#     direction = -1  # 反转颜色顺序以匹配反转后的图例
#   ) +
#   scale_y_continuous(expand = c(0, 0))
ggsave("./figs/图2 分年份各职业人群发病构成比.png", p_各年龄组发病情况及占比)

## 1.3 职业特征------

### 1.3.1 总体职业分布------
occupation_dist <- calculate_distribution(data, c("population_cat"))

print("总体职业分布:")
print(occupation_dist)
writexl::write_xlsx(occupation_dist, "./output/表格1.3.1_人群分布_总体职业.xlsx")

if (!"职业_总体" %in% names(wb)) {addWorksheet(wb, "职业_总体")}
writeData(wb, sheet = "职业_总体", x = occupation_dist)

### 1.3.2 分年份各职业人群发病情况（发病数）及占比------
yearly_occupation_dist <- calculate_distribution(data, c("year", "population_cat")) %>%
  arrange(year, population_cat)
writexl::write_xlsx(yearly_occupation_dist, "./output/表格1.3.2_人群分布_分年份各职业.xlsx")

if (!"职业_逐年" %in% names(wb)) {addWorksheet(wb, "职业_逐年")}
writeData(wb, sheet = "职业_逐年", x = yearly_occupation_dist)


# 绘制图2: 分年份各职业人群发病构成比
yearly_occupation_dist$population_cat %>% unique()

# 1.3.2绘制图2步骤1: 提取2023年的职业发病率排序（排除Others）
# order_2023 <- yearly_occupation_dist %>%
#   filter(year == 2023, population_cat != "Others") %>% 
#   arrange(desc(percentage)) %>% 
#   pull(population_cat)

order_2023 <- c("Farmers" ,"Livestock Related Occupations","Students and Children","Freelancer" ,                  
                "Cadres and Staff" ,"Workers"   , "Retirees"     , "Commercial Services" )


# 1.3.2绘制图2步骤2: 将Others添加到排序末尾
final_order <- c(order_2023, "Others")

# 1.3.2步骤3: 将分类变量转换为因子并固定顺序
yearly_occupation_dist1 <- yearly_occupation_dist %>%
  mutate(population_cat = factor(
    population_cat,
    levels = rev(final_order)))  # 反转顺序以实现大值在底部
  
# 1.3.2绘制图2步骤4: 绘图代码
p_各职业人群发病构成比 <- ggplot(yearly_occupation_dist1, aes(x = factor(year), y = percentage, fill = population_cat)) +
  geom_bar(stat = "identity") +
  labs(
    # title = "各年份各职业人群发病构成比（按2023年比例排序）",
    x = "Year",
    y = "Percentage (%)",
    fill = "Occupation Group"  ) +
  theme_minimal() +
  scale_fill_manual(
    values = setNames(brewer.pal(length(final_order), "Set3"), final_order),  # 精准颜色匹配
    breaks = final_order,  # 图例顺序从高到低
    guide = guide_legend(reverse = F)  ) +  # 图例反转以匹配堆叠顺序
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
    legend.text = element_text(size = 9),
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5, face = "bold") )
ggsave("./figs/图2 分年份各职业人群发病构成比.png", p_各职业人群发病构成比)
# ggsave(p_各职业人群发病构成比,"./figs/图2 分年份各职业人群发病构成比.png")

### 1.3.3 不同职业的性别比------
occupation_gender_ratio <- data %>%
  group_by(population_cat, sex) %>%
  summarise(count = n(), .groups = 'drop') %>%
  group_by(population_cat) %>%
  mutate(total = sum(count),
         ratio = ifelse(sex == "male", count / lag(count), NA)) %>%
  filter(!is.na(ratio))

print("不同职业的性别比:")
print(occupation_gender_ratio)
writexl::write_xlsx(occupation_gender_ratio, "./output/表格1.3.3_人群分布_不同职业性别比.xlsx")

if (!"职业_分性别" %in% names(wb)) {addWorksheet(wb, "职业_分性别")}
writeData(wb, sheet = "职业_分性别", x = occupation_gender_ratio)


### 1.3.4新增 特殊职业分布（牧民及相关工作）------
data_livestock <- data %>% filter(population_cat %in% c("Livestock Related Occupations"))
occupation_dist_livestock <- calculate_distribution(data_livestock, c("pc_cn0"))

print("牧民及相关工作体职业分布:")
print(occupation_dist_livestock) #有待于优化pc_cn0分类
writexl::write_xlsx(occupation_dist_livestock, "./output/表格1.3.4_人群分布_牧民及相关工作.xlsx")

if (!"职业_牧民相关" %in% names(wb)) {addWorksheet(wb, "职业_牧民相关")}
writeData(wb, sheet = "职业_牧民相关", x = occupation_dist_livestock)

## 1.4 确诊病例与临床诊断病例的发病数及构成比变化趋势分析------

case_category_trend <- data %>%
  # mutate(year = year(onset_date)) %>%
  calculate_distribution(c("year","case_category")) %>%
  arrange(year)
writexl::write_xlsx(case_category_trend, "./output/表格1.4_人群分布_确诊病例与临床诊断.xlsx")

if (!"确诊vs临床" %in% names(wb)) {addWorksheet(wb, "确诊vs临床")}
writeData(wb, sheet = "确诊vs临床", x = case_category_trend)

# 构成比
# p_确诊病例与临床诊断病例 <- ggplot(case_category_trend, aes(x = factor(year), y = percentage, fill = case_category)) +
#   geom_bar(stat = "identity") +
#   labs(title = "确诊病例与临床诊断病例的变化趋势",
#        x = "Year",
#        y = "Percentage (%)")
# ggsave("./figs/图1.4 确诊病例与临床诊断病例构成比.png",p_确诊病例与临床诊断病例)
# 更新标签为英文
case_category_trend <- case_category_trend %>%
  mutate(case_category = case_when(
    case_category == "确诊病例" ~ "Confirmed Cases",
    case_category == "临床诊断病例" ~ "Clinically Diagnosed Cases",
    TRUE ~ case_category
  ))

# 绘图

p_确诊病例与临床诊断病例 <- ggplot(case_category_trend, aes(x = factor(year), y = percentage, fill = case_category)) +
  geom_bar(stat = "identity") +
  labs(
    # title = "Trends in the Percentages of Confirmed and Clinically Diagnosed Cases", 
    x = "Year", y = "Percentage (%)", fill = "Case Category") +
  theme_minimal() +
  scale_fill_manual(values = brewer.pal(9, "Set3")[c(1,4)]) +  # 使用Set3色系
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10), 
        legend.text = element_text(size = 9), 
        legend.position = "bottom",
        plot.title = element_text(hjust = 0.5, face = "bold"))

# 保存
if (!dir.exists("./figs")) dir.create("./figs")
ggsave("./figs/Figure_1.4_Trends_in_Cases.png", p_确诊病例与临床诊断病例)

# 发病数
ggplot(case_category_trend, aes(x = factor(year), y = count, fill = case_category)) +
  geom_bar(stat = "identity") +
  labs(
    # title = "Trends in the Count of Confirmed and Clinically Diagnosed Cases",
       x = "Year",y = "Count", fill = "Case Category")+
  theme_minimal() +
  theme(legend.position = "bottom")
  scale_fill_manual(values = brewer.pal(9, "Set3")[c(1,4)])   # 使用Set3色系
ggsave("./figs/图1.4 确诊病例与临床诊断病例的变化趋势.png")

## 1.5 诊断-延误：发病至诊断时间情况（中位数，四分位间距）有问题------
psych::describe(data$diag_delay)
summary(data$onset_date)
summary(data$diag_time)

#下面这些人的发病时间有问题【待讨论】

data %>% filter(onset_date <= as.Date("2003-01-01")) #9人
summary(data$year);summary(data$diag_delay) 
data_debug <- data %>% filter(onset_date >= as.Date("2005-01-01"))
# data_debug <- data %>% filter(diag_delay <=365)
summary(data_debug$diag_delay) #清理后，最大诊断延误时间为3901天
# data <- data %>% dplyr::select(card_id, onset_date:month) 

# 定义一个函数来移除异常值
remove_outliers <- function(df, column) {
  qnt <- quantile(df[[column]], probs=c(.25, .75), na.rm = T)
  iqr <- IQR(df[[column]], na.rm = T)
  df[!(df[[column]] < (qnt[1] - 5*iqr) | df[[column]] > (qnt[2] + 5*iqr)), ]#20250321改
}

# 应用函数到数据框
data_cleaned <- remove_outliers(data, "diag_delay")
summary(data_cleaned$diag_delay)

# 绘制没有异常值的箱式图
# 定义颜色
if (TRUE) {  colors <- c("#8DD3C7", "#FFFFB3", "#BEBADA", "#FB8072", "#80B1D3", "#FDB462", "#B3DE69", "#FCCDE5", "#D9D9D9", "#BC80BD", "#CCEBC5", "#FFED6F")
  print(colors)}

# p_诊断延误 <- ggplot(data_cleaned, aes(x = factor(year), y = diag_delay)) +# data_debug # data_cleaned
#   geom_boxplot() +
#   labs(
#     title = "诊断延迟时间 (diag_delay) 的年度分布 (移除异常值)",
#     x = "年份",
#     y = "诊断延迟天数"
#   ) +
#   theme_minimal() +
#   theme(
#     axis.text.x = element_text(angle = 45, hjust = 1),
#     plot.title = element_text(hjust = 0.5)
#   )

p_诊断延误 <- ggplot(data_cleaned, aes(x = factor(year), y = diag_delay)) +
  geom_boxplot(outlier.shape = NA, fill = "#8DD3C7", color = "#FB8072") +
  # geom_jitter(width = 0.2, alpha = 0.1, color = "#80B1D3") +
  labs(
    # title = "Annual Distribution of Diagnosis Delay (Values over 60 not shown)",
    x = "Year", y = "Diagnosis Delay (Days)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10), 
        legend.position = "bottom",
        plot.title = element_text(hjust = 0.5, face = "bold"), text = element_text(family = "sans")) +
  coord_cartesian(ylim = c(NA, 90))  # 限制 y 轴范围为 [最小值, 50]
ggsave("./figs/图1.5 发病至诊断时间箱式图.png", p_诊断延误)



### 1.5.1 整体发病至诊断时间情况（中位数，四分位间距）------
# 2015年1月26日更新：将data更换为清晰了发病-诊断时间的data_cleaned
overall_diag_time <- data %>% #filter(onset_date >= as.Date("2004-01-01")) %>% 
  summarise(
    median_days = median(as.numeric(difftime(diag_time, onset_date, units = "days")), na.rm = TRUE),
    iqr_days = IQR(as.numeric(difftime(diag_time, onset_date, units = "days")), na.rm = TRUE),
    p25_diag_delay = quantile(as.numeric(difftime(diag_time, onset_date, units = "days")), 0.25, na.rm = TRUE),
    p75_diag_delay = quantile(as.numeric(difftime(diag_time, onset_date, units = "days")), 0.75, na.rm = TRUE)
  )

print("整体发病至诊断时间情况（中位数，四分位间距）:")
print(overall_diag_time)
writexl::write_xlsx(overall_diag_time, "./output/表格1.5.1_发病至诊断时间_整体.xlsx")

if (!"诊断延误" %in% names(wb)) {addWorksheet(wb, "诊断延误")}
writeData(wb, sheet = "诊断延误", x = overall_diag_time)

### 1.5.2 全省逐年发病至诊断时间情况（中位数，四分位间距）------
yearly_diag_time <- data %>%
  group_by(year) %>%
  summarise(
    median_days = median(as.numeric(difftime(diag_time, onset_date, units = "days")), na.rm = TRUE),
    iqr_days = IQR(as.numeric(difftime(diag_time, onset_date, units = "days")), na.rm = TRUE),
    min__diag_delay = min(as.numeric(difftime(diag_time, onset_date, units = "days")), na.rm = TRUE),
    p25_diag_delay = quantile(as.numeric(difftime(diag_time, onset_date, units = "days")), 0.25, na.rm = TRUE),
    p75_diag_delay = quantile(as.numeric(difftime(diag_time, onset_date, units = "days")), 0.75, na.rm = TRUE),
    mx__diag_delay = max(as.numeric(difftime(diag_time, onset_date, units = "days")), na.rm = TRUE),
    .groups = 'drop'
  )

print("逐年发病至诊断时间情况（中位数，四分位间距）:")
print(yearly_diag_time) # 2007年的中卫发病-诊断时间很奇怪！ 2025-01-26;应该是极端值较多？
writexl::write_xlsx(yearly_diag_time, "./output/表格1.5.2_发病至诊断时间_全省逐年.xlsx")

if (!"诊断延误_逐年" %in% names(wb)) {addWorksheet(wb, "诊断延误_逐年")}
writeData(wb, sheet = "诊断延误_逐年", x = yearly_diag_time)

### 1.5.3 地州逐年发病至诊断时间情况（中位数，四分位间距）------
city_yearly_diag_time <- data %>%
  mutate(year = year(onset_date)) %>%
  filter(year >= 2004) %>% #2005年3月21日新增
  group_by(year,city_code) %>%
  summarise(
    median_days = median(as.numeric(difftime(diag_time, onset_date, units = "days")), na.rm = TRUE),
    iqr_days = IQR(as.numeric(difftime(diag_time, onset_date, units = "days")), na.rm = TRUE),
    p25_diag_delay = quantile(as.numeric(difftime(diag_time, onset_date, units = "days")), 0.25, na.rm = TRUE),
    p75_diag_delay = quantile(as.numeric(difftime(diag_time, onset_date, units = "days")), 0.75, na.rm = TRUE),
    .groups = 'drop') %>% arrange(city_code) %>% 
  left_join(nameCE_city)
city_yearly_diag_time$city_code %>% unique();nameCE_city$city_code
city_yearly_diag_time$year %>% unique()

# names(data)
print("地州逐年发病至诊断时间情况（中位数，四分位间距）:")

writexl::write_xlsx(city_yearly_diag_time, "./output/表格1.7_City_yearly_diag_time.xlsx")


if (!"诊断延误_地州" %in% names(wb)) {addWorksheet(wb, "诊断延误_地州")}
writeData(wb, sheet = "诊断延误_地州", x = city_yearly_diag_time)

# 各地州的延迟诊断图
pop <- read_excel("data/origin/发病率2004-2023年v1.xlsx", 
                  sheet = "Report(基于六普七普人口数)", col_types = c("text", "text", "text", 
                                                             "skip", "skip", "numeric", "numeric", "numeric", "numeric", "numeric", 
                                                             "numeric", "numeric", "numeric","numeric", "numeric", "numeric", "numeric", 
                                                             "numeric", "numeric","numeric", "numeric", "numeric", "numeric", "numeric", "skip")) %>%
  filter(region != "county兵团")
name_city <- pop %>% pivot_longer(cols = -c(地区,地区编码,region),names_to = "year") %>%
  filter(region == "city") %>%  ungroup() %>%
  dplyr::mutate(地区编码1=substr(地区编码,1,4)) %>%
  dplyr::select(地区,地区编码=地区编码1) %>% unique()

data1 <- data %>% left_join(name_city,by = c("city_code" = "地区编码")) %>%
  left_join(nameCE_city)


p_诊断延误_分地州 <- ggplot(data1, aes(x = factor(year), y = diag_delay)) +
  geom_boxplot(outlier.shape = NA, fill = "#8DD3C7", color = "#FB8072") +  # 箱式图，填充色和边框色与参考代码一致
  facet_wrap(~ nameE, scales = "free_y") +  # 分面绘制，按地区分面
  labs(
    title = "Diagnosis Delay by City and Year",
    x = "Year",
    y = "Diagnosis Delay (Days)"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),  # 调整x轴标签角度
    plot.title = element_text(hjust = 0.5, face = "bold"),  # 标题居中并加粗
    legend.position = "none",  # 不显示图例
    text = element_text(family = "sans")  # 字体与参考代码一致
  ) +
    coord_cartesian(ylim = c(NA, 300))  # 限制 y 轴范围为 [最小值, 50]
ggsave("./figs/图1.5 发病至诊断时间箱式图（分地州）.png", p_诊断延误_分地州)

### 1.5.4 区县逐年发病至诊断时间情况（中位数，四分位间距）------
county_yearly_diag_time <- data %>%
  mutate(year = year(onset_date)) %>%
  group_by(year,county_code) %>%
  summarise(
    median_days = median(as.numeric(difftime(diag_time, onset_date, units = "days")), na.rm = TRUE),
    iqr_days = IQR(as.numeric(difftime(diag_time, onset_date, units = "days")), na.rm = TRUE),
    p25_diag_delay = quantile(as.numeric(difftime(diag_time, onset_date, units = "days")), 0.25, na.rm = TRUE),
    p75_diag_delay = quantile(as.numeric(difftime(diag_time, onset_date, units = "days")), 0.75, na.rm = TRUE),
    .groups = 'drop'
  )
# names(data)
print("地州逐年发病至诊断时间情况（中位数，四分位间距）:")
print(yearly_diag_time)
writexl::write_xlsx(county_yearly_diag_time, "./output/表格1.8_county_yearly_diag_time.xlsx")


if (!"诊断延误_区县" %in% names(wb)) {addWorksheet(wb, "诊断延误_区县")}
writeData(wb, sheet = "诊断延误_区县", x = county_yearly_diag_time)



# # 检查2004年的异常数据
# data %>% filter(year <= 2004) %>% select(card_id, onset_date, diag_time)

### 1.5.5 区县逐年发病至诊断时间情况（中位数，四分位间距）------

# 0. 
all_summary <- data %>%
  group_by(city_code, county_code, year, month, sex, population_cat) %>%
  summarise(
    min = min(diag_delay, na.rm = TRUE),
    max = max(diag_delay, na.rm = TRUE),
    mean = mean(diag_delay, na.rm = TRUE),
    median = median(diag_delay, na.rm = TRUE),
    p25 = quantile(diag_delay, 0.25, na.rm = TRUE),
    p75 = quantile(diag_delay, 0.75, na.rm = TRUE),
    count_over_180 = sum(diag_delay > 180, na.rm = TRUE),
    .groups = 'drop'
  )
all_summary
writexl::write_xlsx(all_summary, path = "./output/表格1.5.3_发病诊断时间_所有类别.xlsx")


if (!"诊断延误_区县_iqr" %in% names(wb)) {addWorksheet(wb, "诊断延误_区县_iqr")}
writeData(wb, sheet = "诊断延误_区县_iqr", x = all_summary)



# 1. 筛选出发病-报告时间间隔超出180天的个案记录
# data_180 <- data %>%
#   filter(diag_delay > 180)
# write_xlsx(data_180, path = "./output/1.5.5.1 发病诊断时间_超180天.xlsx")

# 2. 统计不同月份（month）的发病-报告时间间隔（diag_delay）的分布情况
month_summary <- data %>%
  group_by(month) %>%
  summarise(
    min = min(diag_delay, na.rm = TRUE),
    max = max(diag_delay, na.rm = TRUE),
    mean = mean(diag_delay, na.rm = TRUE),
    median = median(diag_delay, na.rm = TRUE),
    p25 = quantile(diag_delay, 0.25, na.rm = TRUE),
    p75 = quantile(diag_delay, 0.75, na.rm = TRUE),
    count_over_180 = sum(diag_delay > 180, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  mutate(group = "Month")
month_summary


if (!"诊断延误_月份" %in% names(wb)) {addWorksheet(wb, "诊断延误_月份")}
writeData(wb, sheet = "诊断延误_月份", x = month_summary)

# 3. 统计各年份的发病-报告时间间隔（diag_delay）的分布情况
year_summary <- data %>%
  group_by(year) %>%
  summarise(
    min = min(diag_delay, na.rm = TRUE),
    max = max(diag_delay, na.rm = TRUE),
    mean = mean(diag_delay, na.rm = TRUE),
    median = median(diag_delay, na.rm = TRUE),
    p25 = quantile(diag_delay, 0.25, na.rm = TRUE),
    p75 = quantile(diag_delay, 0.75, na.rm = TRUE),
    count_over_180 = sum(diag_delay > 180, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  mutate(group = "Year")
year_summary # 2007年p25是0???

if (!"诊断延误_年份" %in% names(wb)) {addWorksheet(wb, "诊断延误_年份")}
writeData(wb, sheet = "诊断延误_年份", x = year_summary)

# 4. 统计各城市（city_code）的发病-报告时间间隔（diag_delay）的分布情况
city_summary <- data %>%
  group_by(city_code) %>%
  summarise(
    min = min(diag_delay, na.rm = TRUE),
    max = max(diag_delay, na.rm = TRUE),
    mean = mean(diag_delay, na.rm = TRUE),
    median = median(diag_delay, na.rm = TRUE),
    p25 = quantile(diag_delay, 0.25, na.rm = TRUE),
    p75 = quantile(diag_delay, 0.75, na.rm = TRUE),
    count_over_180 = sum(diag_delay >= 180, na.rm = TRUE),
    count_over_90_180 = sum(diag_delay >= 90 & diag_delay < 180, na.rm = TRUE),
    count_over_15_90 = sum(diag_delay >=15 & diag_delay < 90, na.rm = TRUE),
    count_over_0_15 = sum(diag_delay >=0 & diag_delay < 15, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  mutate(group = "City")
city_summary

if (!"诊断延误_地州_摘要" %in% names(wb)) {addWorksheet(wb, "诊断延误_地州_摘要")}
writeData(wb, sheet = "诊断延误_地州_摘要", x = city_summary)


# 5. 统计不同性别（sex）的发病-报告时间间隔（diag_delay）的分布情况
gender_summary <- data %>%
  group_by(sex) %>%
  summarise(
    min = min(diag_delay, na.rm = TRUE),
    max = max(diag_delay, na.rm = TRUE),
    mean = mean(diag_delay, na.rm = TRUE),
    median = median(diag_delay, na.rm = TRUE),
    p25 = quantile(diag_delay, 0.25, na.rm = TRUE),
    p75 = quantile(diag_delay, 0.75, na.rm = TRUE),
    count_over_180 = sum(diag_delay >= 180, na.rm = TRUE),
    count_over_90_180 = sum(diag_delay >= 90 & diag_delay < 180, na.rm = TRUE),
    count_over_15_90 = sum(diag_delay >=15 & diag_delay < 90, na.rm = TRUE),
    count_over_0_15 = sum(diag_delay >=0 & diag_delay < 15, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  mutate(group = "Gender")
gender_summary

if (!"诊断延误_性别" %in% names(wb)) {addWorksheet(wb, "诊断延误_性别")}
writeData(wb, sheet = "诊断延误_性别", x = gender_summary)

# 6. 统计不同职业分类（population_cat）的发病-报告时间间隔（diag_delay）的分布情况
occupation_summary <- data %>%
  group_by(population_cat) %>%
  summarise(
    min = min(diag_delay, na.rm = TRUE),
    max = max(diag_delay, na.rm = TRUE),
    mean = mean(diag_delay, na.rm = TRUE),
    median = median(diag_delay, na.rm = TRUE),
    p25 = quantile(diag_delay, 0.25, na.rm = TRUE),
    p75 = quantile(diag_delay, 0.75, na.rm = TRUE),
    count_over_180 = sum(diag_delay >= 180, na.rm = TRUE),
    count_over_90_180 = sum(diag_delay >= 90 & diag_delay < 180, na.rm = TRUE),
    count_over_15_90 = sum(diag_delay >=15 & diag_delay < 90, na.rm = TRUE),
    count_over_0_15 = sum(diag_delay >=0 & diag_delay < 15, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  mutate(group = "Occupation")
occupation_summary

if (!"诊断延误_职业" %in% names(wb)) {addWorksheet(wb, "诊断延误_职业")}
writeData(wb, sheet = "诊断延误_职业", x = occupation_summary)

# 7. 统计不同区县（county_code）的发病-报告时间间隔（diag_delay）的分布情况
county_summary <- data %>%
  group_by(county_code) %>%
  summarise(
    min = min(diag_delay, na.rm = TRUE),
    max = max(diag_delay, na.rm = TRUE),
    mean = mean(diag_delay, na.rm = TRUE),
    median = median(diag_delay, na.rm = TRUE),
    p25 = quantile(diag_delay, 0.25, na.rm = TRUE),
    p75 = quantile(diag_delay, 0.75, na.rm = TRUE),
    count_over_180 = sum(diag_delay >= 180, na.rm = TRUE),
    count_over_90_180 = sum(diag_delay >= 90 & diag_delay < 180, na.rm = TRUE),
    count_over_15_90 = sum(diag_delay >=15 & diag_delay < 90, na.rm = TRUE),
    count_over_0_15 = sum(diag_delay >=0 & diag_delay < 15, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  mutate(group = "County")
county_summary


if (!"诊断延误_区县" %in% names(wb)) {addWorksheet(wb, "诊断延误_区县")}
writeData(wb, sheet = "诊断延误_区县", x = county_summary)
saveWorkbook(wb, file = "./output/0_人群分布特征_逐年.xlsx", overwrite = TRUE)

# 创建一个包含所有数据框的列表，每个数据框作为一个元素
data_list <- list(
  month_summary = month_summary,
  year_summary = year_summary,
  city_summary = city_summary,
  gender_summary = gender_summary,
  occupation_summary = occupation_summary,
  county_summary = county_summary
)

# 将列表中的所有数据框导出到同一个 Excel 文件的不同工作表中
write_xlsx(data_list, path = "./output/表格1.5.4_发病诊断时间_分类别.xlsx")

# 2. 人群分布特征分析（分四阶段）------
head(data)
names(data)
# 定义年份阶段
year_periods <- list(
  "All Years" = range(data$year),
  "2005-2009" = c(2005, 2009),
  "2010-2015" = c(2010, 2015),
  "2016-2019" = c(2016, 2019),
  "2020-2023" = c(2020, 2023)
)

# 函数：根据给定的年份范围进行分组并计算统计数据
compute_stats <- function(data, year_range) {
  filtered_data <- data %>%
    filter(year >= year_range[1] & year <= year_range[2])
  
  total_cases <- nrow(filtered_data)
  
  age_group_cases <- filtered_data %>%
    count(age_group, name = "count") %>%
    mutate(category = "Age Group", subgroup = age_group) %>%
    select(category, subgroup, count)
  
  sex_cases <- filtered_data %>%
    count(sex, name = "count") %>%
    mutate(category = "Sex", subgroup = sex) %>%
    select(category, subgroup, count)
  
  population_cat_cases <- filtered_data %>%
    count(population_cat, name = "count") %>%
    mutate(category = "Population Category", subgroup = population_cat) %>%
    select(category, subgroup, count)
  
  patient_sources_cases <- filtered_data %>%
    count(patient_sources, name = "count") %>%
    mutate(category = "Patient Sources", subgroup = patient_sources) %>%
    select(category, subgroup, count)
  # 
  # city_code_cases <- filtered_data %>%
  #   count(city_code, name = "count") %>%
  #   mutate(category = "City Code", subgroup = city_code) %>%
  #   select(category, subgroup, count)
  
  # 总病例数作为一个单独的行
  total_df <- tibble(category = "Total", 
                     subgroup = "Total", 
                     count = total_cases)
  
  # 合并所有统计数据
  bind_rows(total_df, age_group_cases, sex_cases, population_cat_cases, patient_sources_cases) %>% #, city_code_cases
    mutate(year_period = paste(year_range[1], "-", year_range[2]))
}

# 应用函数到每个年份阶段，并存储结果
results_list <- lapply(names(year_periods), function(period_name) {
  period_range <- year_periods[[period_name]]
  compute_stats(data, period_range)
})

# 将所有结果合并成一个数据框
results <- bind_rows(results_list)

# 将数据转为宽格式以便展示
results_wide <- results %>%
  pivot_wider(names_from = year_period, values_from = count, names_prefix = "Y")

# 打印结果
print(results_wide)

# 如果需要更美观的输出，可以使用kable
# library(knitr)
# kable(results_wide, format = "markdown")

write_xlsx(results_wide, path = "./output/Table1_人群分布特征分析（分阶段）.xlsx")

# # 3. 人群分布特征分析（原先的分三阶段；0227）------
# head(data)
# names(data)
# # 定义年份阶段
# year_periods <- list(
#   "All Years" = range(data$year),
#   "2005-2015" = c(2005, 2015),
#   "2016-2020" = c(2016, 2020),
#   "2021-2023" = c(2021, 2023)
# )
# 
# # 函数：根据给定的年份范围进行分组并计算统计数据
# compute_stats <- function(data, year_range) {
#   filtered_data <- data %>%
#     filter(year >= year_range[1] & year <= year_range[2])
#   
#   total_cases <- nrow(filtered_data)
#   
#   age_group_cases <- filtered_data %>%
#     count(age_group, name = "count") %>%
#     mutate(category = "Age Group", subgroup = age_group) %>%
#     select(category, subgroup, count)
#   
#   sex_cases <- filtered_data %>%
#     count(sex, name = "count") %>%
#     mutate(category = "Sex", subgroup = sex) %>%
#     select(category, subgroup, count)
#   
#   population_cat_cases <- filtered_data %>%
#     count(population_cat, name = "count") %>%
#     mutate(category = "Population Category", subgroup = population_cat) %>%
#     select(category, subgroup, count)
#   
#   patient_sources_cases <- filtered_data %>%
#     count(patient_sources, name = "count") %>%
#     mutate(category = "Patient Sources", subgroup = patient_sources) %>%
#     select(category, subgroup, count)
#   # 
#   # city_code_cases <- filtered_data %>%
#   #   count(city_code, name = "count") %>%
#   #   mutate(category = "City Code", subgroup = city_code) %>%
#   #   select(category, subgroup, count)
#   
#   # 总病例数作为一个单独的行
#   total_df <- tibble(category = "Total", 
#                      subgroup = "Total", 
#                      count = total_cases)
#   
#   # 合并所有统计数据
#   bind_rows(total_df, age_group_cases, sex_cases, population_cat_cases, patient_sources_cases) %>% #, city_code_cases
#     mutate(year_period = paste(year_range[1], "-", year_range[2]))
# }
# 
# # 应用函数到每个年份阶段，并存储结果
# results_list <- lapply(names(year_periods), function(period_name) {
#   period_range <- year_periods[[period_name]]
#   compute_stats(data, period_range)
# })
# 
# # 将所有结果合并成一个数据框
# results <- bind_rows(results_list)
# 
# # 将数据转为宽格式以便展示
# results_wide <- results %>%
#   pivot_wider(names_from = year_period, values_from = count, names_prefix = "Y")
# 
# # 打印结果
# print(results_wide)
# 
# # 如果需要更美观的输出，可以使用kable
# # library(knitr)
# # kable(results_wide, format = "markdown")
# 
# write_xlsx(results_wide, path = "./output/Table1_人群分布特征分析（分阶段）.xlsx")
